1 Background

This is a tutorial on how to process your 16S ITS1/ITS2 amplicon sequences and identify the taxonomic identification of the ASVs (i.e., amplicon sequence variants, also known as zOTUs for zero OTUs or ESVs for exact sequence variants) in your sequence data.

To create this tutorial, I have assembled scripts I’ve used to analyze 16S data provided by Dr. Emma Aronson’s lab. The data I am working with to create this workflow comes from a project that examined soil microbial community composition in Mount Saint Helens. The target region was the V4 region within the 16S gene, and sequencing was performed with an Illumina MiSeq (2x300).

This tutorial would not have been possible without Dr. Benjamin Callahan’s DADA2 program (Callahan et al. 2016) and tutorials. Additionally, I would like to especially thank Dr. Mike Lee for his guidance, his patience, and his Happy Belly Bioinformatics tutorial called Amplicon Analysis tutorial (Lee 2019).

1.1 Considerations before you begin

I was able to analyze these sequences on a High Performance Computing cluster (HPCC) that uses a Slurm scheduler. The minimum amount of total memory I used (not per CPU, but overall) for each step in this workflow (i.e., each step as a separate ‘job’ via Slurm) was 400GB. Having enough memory is essential to running most of these programs, so please keep this in mind before embarking on this workflow!

These steps are also time consuming depending on your memory constraints, so do not be concerned if this process takes a while. If you plan to run through multiple steps in the workflow in one sitting, then I suggest loading tmux before you run your scripts. Here is a handy tmux cheat sheet that I refer to often. For more information on what tmux is and how to utilize it, check this tmux crash course by Josh Clayton.

I also suggest exploring a program called neovim (aka nvim) that allows you to use Vim (a text editor) to edit R code and run the code simultaneously. Though nvim is not necessary to run through this workflow, I find that it makes my life a bit easier when running through the DADA2 portion of the workflow. I will get more into the usage of nvim once we get to the DADA2 step(s), but for more information please view the Neovim Github as well as its documentation. You can also find a helpful nvim tutorial here created by Dr. Thomas Girke from UC Riverside.

Additionally, you will need to change your path to each of these programs depending on where they are stored in your computer or HPCC. If you are running these steps locally (which, if you are, then you have one badass computer!), then you can skip the module loading lines in each step – loading modules is specifically for running these scripts on a HPCC that uses a Slurm Workload Manager.

1.2 Submitting Scripts as Jobs with Slurm

If you are unsure as to how to set up the script for submitting on your HPCC, check the code chunk below. This is the information I use when submitting a job to our Slurm system. Again, this is specifically for a system that uses the Slurm scheduler. For more information on what the arguments mean and how to set up your job submissions, please refer to this handy cheatsheet made by Slurm.

When I don’t know exactly what a program’s output will look like, I will run the program via an interactive job on the HPCC. I also suggest running programs interactively if the program requires multiple lines of code to run, and you want to make sure each step has the correct input (whether it be a file, an object, or the output of a previous step in the code). For some more information on interactive jobs in a Slurm system, check out this blog post by Yunming Zhang. This is how I set up an interactive job on the HPCC (that uses Slurm).

1.3 A bash scripting tip for before we start

I wanted to share a bit of code that you will see being implemented in every script throughout the tutorial. This little bit of code will help you pull out the sample names from your files, allowing you easily run through your files while also keeping track of which samples those files belong to.

Here I am using using a for loop to loop through each fastq files in a specific directory. In each iteration of the loop, an $f variable is created, which uses the basename function to get the file name of the $FILE variable. Then $SAMPLE is created by using % to remove the .fastq extension and everything that follows, keeping only the file name (minus the extension) and calling that $SAMPLE. Then we can use the $SAMPLE variable to substitute the file names, which come in handy for running these scripts over multiple samples at one time. This concept should become clearer as we move through the workflow. If you’d like more information on string substitution (i.e., using % to remove parts of a string), please see this helpful link.

2 Sample pre-processing

2.1 Demultiplex your samples

When preparing sequencing libraries, we typically multiplex our samples. This means that during library preparation, we’ve attached barcodes to our sequences that help us trace the sample that these sequences came from. This allows us to pool multiple libraries together in one sequencing run. After sequencing, the sequences are demultiplexed, meaning the individual sequences are separated out by sample into individual FASTQ files.

Typically your samples will be returned to you already demultiplexed. However, if your samples are still pooled into one large FASTQ file, do not panic! You can follow the demultiplexing tutorial by Dr. Mike Lee which utilizes the sabre tool. Or, you can use bcl2fastq2 by Illumina (more information here).

2.2 Sequence Quality and Where to Trim

2.2.1 Check the quality of your sequences with FastQC

It’s always a good idea to check the quality of your sequences before you start your analysis, regardless of the type of sequences they are (metagenomes, RNA-seq data, etc). FastQC (Andrews, n.d.) provides a comprehensive report on the quality of your sequences and is helpful for the following: identifying primers or adapters still attached to your sequences; determining the quality of your reverse reads; etc. You can also use the FastQC reports to determine if you should attempt to merge your forward and reverse reads, or just proceed with only the forward reads.

FastQC will return a report assessing the per base and per sequence quality of your sequences, as well as the GC and N (i.e., unidentified base) content across your sequences, the distribution of your sequence lengths, and whether or not adapters are still attached to your sequences. The second tab of the report details the per base sequence quality across all of your sequences. The per base quality score (Q score), also known as a Phred score, is the estimated probability that the base call is wrong. The following equation is used for calculating the Q score: \[ Q = -10log_{10}E \] Here, E is the estimated probability of the base call being wrong. The higher the Q score, the smaller the probability of a base call error. A quality score of 30 (Q30) means that the probability of an incorrect base call is 1 in 1000, and that the base call accuracy (1 - probability of incorrect base call) is 99.9%. For more information on quality scores, please see this info from Illumina.

Below is an example of the “per base sequence quality” portion of the report. This portion of the report helps me to determine where I should trim my sequences as I move forward with the analysis. This part of the report can also give you a sense on whether there was an error in your sequencing run. For example, if the average quality score (i.e., the blue line in the report) across all of the bases dips below 30 for half of the sequence length in all of my samples, that could indicate that there was an error with the sequencing run itself.

Figure 1: Per Base Quality Scores from FastQC Report


Another useful piece of the FastQC report is the adapter content tab, which is the very last tab in the report. This portion of the report tells us the percentage of reads that have adapter sequences at specific base positions along the reads. The following snapshot from a FastQC report shows that the Small RNA 3’ adapter sequence is found in ~2% of the sequences starting at around the 160th base. We can use this information to then decide exactly which adapter sequences to cut from our samples in the trimming step.

Figure 2: Frequency of Adapter Sequences from FastQC Report


For more on how to interpret FastQC reports, please check out this helpful FastQC tutorial from Michigan State University.

2.2.2 Expected Error Filtering of Sequences with eestats

The eestats2 program (Edgar and Flyvbjerg 2015) creates a report detailing the percentage of reads that will pass through an expected error filter when the reads are at different lengths. Specifically the program will determine how many reads at each specific length (i.e., 50 bp, 100 bp, 150 bp, etc.) have good enough quality to surpass the three expected error thresholds: 0.5%, 1%, and 2%.

Before you run the eestats program, be sure to gunzip (aka decompress) your fastq.gz files! You can do that by running the following command: gunzip /path/to/*.fastq.gz.

The following image shows a sample eestats2 report. For this specific sample, the ideal trimming length of the sequences would be around 250 basepairs long. This is because when considering the expected error threshold of 1%, more than 80.8% of the sequence pass this threshold. Though the 300 bp length also allows high rentention of reads, we know that the per base quality of this sample drops as we approach the 300 bp position. Thus, it seems like trimming these sequences to 250 bps would be ideal moving forward.
Figure 3: eestats Report Example


For more information on the eestats2 programs by USEARCH, please read the documentation here.

2.3 Trim sequences and cut adapters, primers, etc

There are plenty of programs out there that can be used for trimming, and the following three are the most popular for amplicon analyses: cutadapt, trimmomatic, and bbduk (Bushnell, n.d.). All of these programs are reputable, but I personally like the bbduk, and will use this tool for trimming and adapter removal.

Before I trim my sequences, I refer to the FastQC reports to find out exactly which adapters I should remove from my sequences. For example, when looking at the adapter content portion of the FastQC report above, I can see that the Nextera Transposase Sequence is still present in that particular sample. Thankfully Illumina shares their adapter sequences on their website, allowing us to easily find common adapters in sequences, like the Nextera Transposase Sequence for example.

I also know that with the sequences I am analyzing, the PCR primers are still attached (FastQC may identify these primers in your report’s Overrepresented Sequences tab, but not necessarily the origin of these sequences). I can either remove these primer sequences using the actual sequence using the (literal=) flag, or I can trim from the right (ftr=) and/or the left (ftl=) of the sequences if I know exactly how long the primer sequences were.

It is recommended to check the overrepresented sequences from the FastQC report to see if there are contaminating sequences present in your data. I suggest taking the most frequent overrepresented sequence and running it through BLASTn if the source of this overrepresented sequence says “No Hit” (meaning that FastQC cannot attribute this sequence to its list of adapter sequences). If the sequence comes up as a contaminant (i.e., a different gene than the amplicon you’re looking at) or adapter/primer of some kind, you can add this to the literal= flag in bbduk to remove the contaminant.

In addition to removing adapter and primer sequences using the the literal= flag, I also include a reference file provided by bbduk (referenced in the ref= flag) that contains all of the Illumina TruSeq adapters. The sequences in the reference file, in addition to the given adapters and primers, will be removed from the sequences. Below the shell script is a description of all of the flags used by bbduk and exactly what they mean.

For more information on the bbduk flags, please see the bbduk documentation.

path=/path/to/sequences/here # replace with the path to your files

if [[ ! -d ./Trimmed_Seqs ]]; then # creating directory to store trimmed sequences in
    mkdir Trimmed_Seqs
fi

if [[ ! -d ./Trimmed_Seqs/16S_Trimmed ]]; then # creating directory for specifically trimmed 16S sequences
    mkdir Trimmed_Seqs/16S_Trimmed
fi

for i in 16S_Seqs/*_R1.fastq;
do
    f=$(basename $i)
    SAMPLE=${f%_R*}
    
    bbduk.sh -Xmx10g in1=${path}/16S_Seqs/${SAMPLE}_R1.fastq in2=${path}/16S_Seqs/${SAMPLE}_R2.fastq out1=${path}/Trimmed_Seqs/16S_Trimmed/${SAMPLE}_R1_clean.fastq out2=${path}/Trimmed_Seqs/16S_Trimmed/${SAMPLE}_R2_clean.fastq literal=TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG,GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG ref=/bigdata/aronsonlab/shared/bbmap_resources/adapters.fa rcomp=t ktrim=r k=23 maq=10 minlength=200 mink=11 hdist=1 tpe tbo
    
done

# ref ---> file provided by bbduk that holds collection of Illumina TruSeq adapters
# literal=(sequence here) ---> literal adapter sequences to remove; "N" represents any base -- in this case, they are indexes within the adapters
# rcomp=t ---> Rcomp looks for kmers and their reverse-complements, rather than just forward kmer, if set to true
# ktrim=r ---> “ktrim=r” is for right-trimming (3′ adapters)
# k=23 ---> look for kmer that is 23 bp long
# mink=11 ---> in addition to kmers of x length, look for shorter kmers with lengths 23 to 11 (in this case)
# maq=10 ---> This will discard reads with average quality below 10
# hdist=1 ---> hamming distance of 1
# mlf=50 ---> (minlengthfraction=50) would discard reads under 50% of their original length after trimming
# trimq=10 ---> quality-trim to Q10 using the Phred algorithm, which is more accurate than naive trimming.
# qtrim=r ---> means it will quality trim the right side only
# tpe ---> which specifies to trim both reads to the same length
# tbo ---> which specifies to also trim adapters based on pair overlap detection using BBMerge (which does not require known adapter sequences)
# mm ----> Maskmiddle ignores the middle base of a kmer, can turn off with mm=f

To be extra cautious and ensure that the trimming step was successful, I will run the trimmed sequences through FastQC and compare the reports. If the per base and per sequence qualities have improved and/or the adapters are absent, then I will move forward with the workflow. However, if I am still not happy with the quality of the trimmed reads, I will then run the trimmed reads through bbduk in hopes of removing persistent, unwanted sequences. I will also check the overrepresented sequences and their frequencies again, and run the most frequent overepresented sequence(s) in BLASTn.

3 ASV Assignment with DADA2

All of the steps in this portion of the workflow (excluding the tmux and nvim code chunks) have been adapted from Dr. Callahan’s DADA2 tutorial and Dr. Lee’s amplicon tutorial.

To prepare for running DADA2, I want to separate our sequence files by locus and region. For example, you do not want to analyze your 16S and ITS2 sequences together in DADA2 – combining loci and even different regions of the same loci can interfere with the DADA2 algorithm. For example, even if you have 16S sequences of just the V3 region, and a set of 16S sequences with the V3-V4 region, you would want to run these regions separately through the DADA2 pipeline. The reason for this will become clearer as we get to the filtering and trimming step and the error rate prediction step.

To create separate directories for your sequence data, I first ensure that their file names include their amplicon that’s been sequenced for that particular sample (e.g., the 16S V4 data for Sample1 is in the file Sample1_16S.V4_R1_001.fastq). Then I would run the following line of code.

Here I am using a * which is a special character that can be used to represent any character or set of characters. In this case, I am telling the mv command to move any files that have the _16S.V4_ pattern anywhere in the file name to a directory called 16S.V4_Seqs. After running this command, I make sure that my script containing the following DADA2 R code is in the directory with the specific files you want to analyze. My DADA2 R script is called DADA2_tutorial_16S_pipeline.R, which you will see me reference in a couple of code chunks.

3.1 Run Interactive Job + tmux on HPCC

Personally, I like to run through DADA2 via an interactive job on our HPCC. This will allow us to run scripts line by line and check the output, rather than submitting a job to run in the cluster without our supervision. Basically, this is an easy to way constantly check our progress and (ideally) catch errors as soon as they happen. Again, your HPCC must use Slurm to run an interactive job in this manner.

Once the interactive job is running, we can use tmux and nvim to start running through the DADA2 R script.

Figure 4a: Neovim Display of R script


Almost immediately after the R script opens in nvim, I type \rf. This will open another window showing your terminal. You can toggle the horizontal verses vertical alignment fo the windows by typing Ctrl-w shift-H for a horizontal alignment or Ctrl-w shift-Vfor a vertical alignment. Below is what the screen should look like after typing \rf followed by Ctrl-w shift-H. You can see the R script is open in the left window, and my terminal is open in the right window.
Figure 4b: View of Neovim Display after typing \rf


Now that we are in nvim, all you need to do to run a line of code is to just hit the space bar! You can also toggle between windows using Ctrl-w w, edit or type code by pressing i to insert code, and leaving the editing mode by pressing esc. To quit and save changes to your R file, just type :wq, or to quit without saving changes to your file, just type :q!. When using nvim, I keep Dr. Girke’s handy tmux/nvim tutorial open as a reference just in case.

3.1.1 Load the path & FASTQ files

We can start by loading the libraries we need as well as the path to the sequences you want to analyze. In this example I will be analyzing 16S V3-V4 sequences, so I set the path object to be the path to those specific sequences.

If you have already started the DADA2 workflow and want to pick up from where you left off, then you can run this next code chunk to load everything that was in your global environment and saved to an .Rdata object (I will explain this code when we save our first .Rdata file).

3.2 Check sequence quality

Though we have already done this with FastQC, there is a step here to check the quality of our forward and reverse reads. These per base sequence quality reports do look nicer than the FastQC output and tell you the total number of reads in that particular sample. We use ggsave() from the ggpubr package to save these reports as high quality PDFs.

Figure 5a: Per Base Sequence Quality Reports for 2 Samples (R1 Only)

Figure 5b: Per Base Sequence Quality Reports for 2 Samples (R2 Only)


3.3 Filter and Trim

Now that we’ve set up our file paths and checked the quality of our sequences, we can set them up for the DADA2 filter and trim step. First, we create objects that will hold the file names of filtered sequences based on the sample names we have provided. Then the filterAndTrim command will filter the reads based upon the following: read quality, read length, the number of Ns (i.e., unknown bases) in a read, the maximum number of expected errors after truncating the reads, and whether or not reads in your sample match the PhiX genome (i.e., a small virus genome used as a control in Illumina sequencing runs; more information here). The maximum expected errors is calculated by solving for E in the Quality Score equation (see above). We have also specified here that we would like our output FASTQ files to be compressed, and that we can multithread this filter and trimming process. Keep in mind that if you are using a Windows to run your analyses locally, then multithreading is not available for this step.

Soem crucial things to consider at this step are your read lengths (i.e., 2x250, 2x300), the locus and region(s) you’re sequencing, the per base quality of your sequences, the and whether you are using paired-end reads or just forward reads.

For the ITS1 and ITS2 genes, their lengths vary so wildly that truncating the sequences based on a specific length is not recommended. However, for the 16S gene, its regions (i.e., V1-V9) are a bit more reliable in their average lengths. The 16S V4 region varies between ~250 - 283 nucleotides in length (Illumina, n.d.; Vargas-Albores et al. 2017), whereas the V3 region varies between ~130 to 190 nucleotides (Vargas-Albores et al. 2017). The V3-V4 region ranges between ~400 - 470 nucleotides in length. For more information, check out this DADA2 GitHub issue as well as Rausch et al. (2019) and Bukin et al. (2019).

Forward and reverse reads are merged if they have at least a 12 base pair overlap. If you are using paired-end reads, then your merged read lengths (considering the 12 nucleotide overlap) need to total up to these region lengths. For example, let’s say you’re truncating your 16S V3-V4 forward reads to 250 base pairs (bp) long and your reverse reads to 160 bp long. If your reads are merged, the total length will be 250 + (160-12) = 398 bp long. This total read length of 398 bp would be a decent minimum read length considering that the range of the 16S V3-V4 region is ~400 - 470 bp.

Lastly, when setting your expected errors per forward and reverse reads (maxEE=c(R1,R2)), it is important to consider the per base sequence quality of your reads. Because reverse reads typically have lower per base sequence quality than your forward reads, you may want to relax the expected errors for your reverse reads.

If few too reads are surviving this step, then consider changing your expected errors per read parameter or adjusting the truncLength of your reads. Referring to your FastQC and eestats2 reports may be provide even more clarity for how you want to define these paramters. For more information on the filterAndTrim function, please view this documentation.

3.3.1 Learn the Error Rates

Dr. Callahan developed an algorithm for a parametric error model that can use both inference and estimations to determine the error rates for each sample. Here is the excerpt of Dr. Callahan describing this function in his DADA2 tutorial.

The DADA2 algorithm makes use of a parametric error model (err) and every amplicon dataset has a different set of error rates. The learnErrors method learns this error model from the data, by alternating estimation of the error rates and inference of sample composition until they converge on a jointly consistent solution. As in many machine-learning problems, the algorithm must begin with an initial guess, for which the maximum possible error rates in this data are used (the error rates if only the most abundant sequence is correct and all the rest are errors).

Figure 6: Rarefaction Curve Example


Once we have constructed the error model from the reads, we can plot the observed frequency of transitions from base to base (i.e., A2A indicates an A followed by an A) as a function of the consensus quality score at that position in the read. The individual points in black represent the observed error rates for each consensus quality score. The black line shows estimated error rates after convergence of the ML algorithm, and the red line shows error rates expected under the nominal definition of the Q-score. One thing we notice is that as consensus quality score increases, the error rates (black lines) decrease, which is expected.

If it seems like we are doing a lot of sanity checks throughout this workflow, it’s because we are! This process can take a while and require some trouble shooting, so it’s good to constantly check your work as you make your way through the workflow. Coding without sanity checks is never recommended.

3.3.2 Identify ASVs in Reads

Before we run the Divisive Amplicon Denoising Algorithm (dada()), we have to remove the files from our filtFs/filtRs objects that were dropped in the filtering step. The algorithm will not work if we include file names of files that do not actually exist in the Filtered directory.

Now it’s time to use the dada algorithm to infer our amplicon sequence variants (ASVs) from our sequences. For more information on how the Divisive Amplicon Denoising Algorithm works, please see Callahan et al. (2016). To increase the signal of sequence variants with very low abundance across samples, you can choose to pool (pool=TRUE) or pseudo-pool (pool=pseudo) your sample sequences together.

The output of the algorithm will be a data-class object, containing the number of ASVs inferred out of the number of total input sequences.

The wonderful thing about ASVs is that because they are assigned based on 99% sequence identity, they are true representative of biological sequences and thus directly comparable across workflows (Prodan et al. 2020; Callahan, McMurdie, and Holmes 2017). I highly recommend reading these papers for more information on the benefits of using ASVs/ESVs/zOTUs.

3.4 Merge Forward + Reverse Reads

At this point in the workflow, we are now going to merge our denoised Forward and Reverse reads to get our contiguous sequences (i.e., contigs). Sequences will be merged if they share at least 12 nucleotides. These sequences must be identical to each other in these overlapping regions or else they will not be merged.

The mergers object is a data.frame containing the merged sequence, its abundance, and some statistics about the sequences themselves. If most of your reads do not merge, then you should revisit the filter and trimming step. It could be that you cut too much off of your sequencing reads, making it more difficult to successfully merge your reads.

3.6 Remove Chimeras

Now we are going to remove all chimeric sequences from our data. Chimeras are the the result of two or more biological sequences incorrectly joining together. This is often a result of PCR for a number of reasons. In DADA2 specifically, chimeras are identified if they can be reconstructed from right and left segments from two or more “parent” sequences. The object seqtab.nochim will be our sequence table with the chimeras removed. Most of your reads should not be removed during this step. However, according to the DADA2 FAQ page if you have more than 25% of your reads removed, then it is likely that primers are still attached to your sequences. Be sure to remove these primers in the trimming step with either bbduk, cutadapt, or trimmomatic and begin the workflow again.

3.7 Track the Reads - Sanity Check

Time for a sanity check and see how many reads we have at this point in our workflow. This is a great place to see if we have lost any reads, and at which steps they were lost - which can really help us determine if we trimmed our reads to the appropriate length. If a lot of reads are lost, it is recommended to check if primers and adapters are still attached to your sequences, and the truncation length of your sequences for the filter and trimming step.

3.8 Assign Taxonomy to ASVs

Now it’s time to assign taxonomic identities to our ASVs. Dr. Callahan utilizes the naive Bayesian classifier his assignTaxonomy function, which is the same classifier used by the Ribosomal Database Project (RDP) for taxonomic assignment. For more information on how this classifier works, please read Wang et al. (2007).

In order to assign the taxonomic IDs to our ASVs, we need to have a reference database FASTA file to use as our known training data. These training data with known references will help the classifier determine which taxa our ASVs belong to. Currently the options for reference databases include the latest versions of the Silva database (for 16S), the Ribosomal Database Project database (for 16S), and the UNITE database (which should be used specifically for ITS sequences). A Green Genes database file is also included (for 16S analyses), but Green Genes has not been updated in a long time and thus is not the best choice for our reference training dataset. Dr. Callahan has included these reference database files here.

If you are only seeing taxonomic identification at the Phyla level (i.e., the rest of the columns are filled with NAs), then this could indicate that we have not trimmed our sequences correctly. For example, let’s say we have reads that are 300bp long (from 2x300 PE sequencing), but we are interested in the 16S V3 region which ranges from ~ 130-190 nucleotides in length. If we have not trimmed our sequences down to our desired region length (here 130-190 nucleotides), then our merged reads are no longer reliable, and the classifier will incorrectly identify our ASVs. Remember, the reads need to have at least a 12 nucleotide overlap to merge - so if we are not trimming our reads correctly, we could create merged sequences that are not accurate representations of the regions we are trying to identify, which will hurt us in the taxonomic assignment step.

3.9 Save DADA2 Output for future analysis

We have finished the DADA2 portion of the workflow! We can save the output from DADA2 as R objects, text files, and tsv files for future import into R.

First we create a vector of the ASV labels called asv_headers that we will use to make the ASV IDs easier to read. We then use a for loop to add an “ASV” prefix to our asv_headers so that they are easily identifiable by an ASV number instead of just a number to represent each ASV. We then combine asv_headers with the sequences themselves to make an object called asv_fasta, which now holds our ASV sequences and their respective IDs.

Now we can save our ASV count table, our ASV taxonomy table, and the ASV sequences themselves as separate files and R objects. I wanted to provide several file options because some people have a preference as to how they import data into R.

# making and writing out a fasta file of our final ASV seqs w/ ASV IDs:
write(asv_fasta, "16S_ASVs_dada2.fa") # write fasta file
write.table(asv_fasta,"16S_ASVs_dada2.txt",sep="\t",row.names=TRUE,col.names=TRUE)

# ASV count table:
asv_counts <- t(seqtab.nochim)
row.names(asv_counts) <- sub(">", "", asv_headers)

# For Vegan format: sample IDs as rows, ASVs as columns
asv_tab<-t(asv_counts)
write.table(asv_tab, "16S.V4_ASVs_Table_dada2.tsv", sep="\t", quote=F, col.names=NA)
write.table(asv_tab,"16S.V4_ASVs_Table_dada2.txt",sep="\t",row.names=TRUE,col.names=TRUE)
# For Phyloseq format: ASVs as row IDs, sample IDs as columns
write.table(asv_counts, "16S.V4_ASVs_Counts_dada2.tsv", sep="\t", quote=F, col.names=NA)
write.table(asv_counts,"16S.V4_ASVs_Counts_dada2.txt",sep="\t",row.names=TRUE,col.names=TRUE)

# taxa ID table:
asv_tax <- taxa
row.names(asv_tax) <- sub(">", "", asv_headers)
write.table(asv_tax, "16S.V4_ASVs_Taxonomy_dada2.tsv", sep="\t", quote=F, col.names=NA)
write.table(asv_tax,"16S.V4_ASVs_Taxonomy_dada2.txt",sep="\t",row.names=TRUE,col.names=TRUE)

#### Save all ASV objects as R objects ####
saveRDS(asv_tax, file = "16S.V4_ASVs_Taxonomy_dada2_Robject.rds", ascii = FALSE, version = NULL,
        compress = TRUE, refhook = NULL)
saveRDS(asv_tab, file = "16S.V4_ASVs_Counts_dada2_Robject.rds", ascii = FALSE, version = NULL,
        compress = TRUE, refhook = NULL)
saveRDS(asv_fasta, file = "16S.V4_ASV_Sequences_dada2_Robject.rds", ascii = FALSE, version = NULL,
        compress = TRUE, refhook = NULL)

Personally, I like using R objects (file extension .rds) in my analyses. In order to import R objects into R, you can run data.frame(readRDS("path/to/Robject.rds", refhook = NULL)) to create a data frame holding the contents of your R object file.

4 Statistical Analysis

At this point you should have either/or R objects, text files, and tsv files containing the following: 1. your ASV sequences in FASTA format, 2. your ASV count table, and 3. your ASV taxonomy table. You should also have some metadata for your samples that will allow for deeper investigation into your microbial data.

4.1 Import and Prepare Data for Analyses

First before we import any data, let’s make sure that we are in the right directory (where our DADA2 files are stored) and that have all of the necessary R libraries loaded.

Now let’s import the DADA2 output into R for some statistical analyses. We will import our ASV count table, our ASV taxonomic table, and our metadata for this dataset. We are also going to create an object called colorset1 to contain the color labels for each of our categories. This will help us keep the colors consistent for each category in all of our figures.

## Import bacterial ASV count data
bac.ASV_counts<-data.frame(readRDS("16S.V4_MSH_ASVs_Counts_dada2_9.20.2021_Robject.rds", refhook = NULL))
dim(bac.ASV_counts)
head(bac.ASV_counts)
colnames(bac.ASV_counts)<-gsub("X1", "1", colnames(bac.ASV_counts)) # shorten sample names to match sample names in metadata file
bac.ASV_counts$ASV_ID<-rownames(bac.ASV_counts)
head(bac.ASV_counts)

## Import metadata
metadata<-as.data.frame(read_excel("MSH_MappingFile_for_Workflow.xlsx"), header=TRUE)
head(metadata)
metadata<-na.omit(metadata) # drop NAs from metadata
head(metadata)
metadata$SampleID<-gsub("(\\_.*?)\\_([0-9])","\\1.\\2", metadata$SampleID) # replace second _ with .
rownames(metadata)<-metadata$SampleID

# create color variable(s) to identify variables by colors
## color for Category
colorset1 = melt(c(ClearCutSoil="#D00000",Gopher="#f8961e",NoGopher="#4ea8de",OldGrowth="#283618"))

colorset1$Category<-rownames(colorset1)
colnames(colorset1)[which(names(colorset1) == "value")] <- "Category_col"
colorset1

metadata<-merge(metadata, colorset1, by="Category")
head(metadata)
metadata$Category_col <- as.character(metadata$Category_col)
rownames(metadata)<-metadata$SampleID

## Import ASV taxonomic data
bac.ASV_tax<-data.frame(readRDS("16S.V4_MSH_ASVs_Taxonomy_dada2_9.20.2021_Robject.rds", refhook = NULL))
head(bac.ASV_tax)

bac.ASV_tax[is.na(bac.ASV_tax)]<- "Unknown" # turn all NAs into "Unkowns"
bac.ASV_tax$Species<-gsub("Unknown", "unknown", bac.ASV_tax$Species) # change uppercase Unkonwn to lowercase unknown for unknown species classification
head(bac.ASV_tax)
bac.ASV_tax$ASV_ID<-rownames(bac.ASV_tax) # create ASV ID column to use for merging data frames
head(bac.ASV_tax)

4.1.1 Data Formatting, Filtering, and Transformation

Now that we have imported our ASV count table, our ASV taxonomy table, and our metadata, we can start to reformat the actual data objects in R to get them ready for running through the vegan suite of tools. First we are going to merge our ASV count table and our ASV taxonomy tables together and filter out some unwanted taxa.

Even though we are analyzing bacterial data, sometimes chloroplast and mitochondrial sequences are mistakenly attributed to 16S genes. For example, in the Silva database, Chloroplast sequences attributed to Eukaryotes are found within the databases’ set of Cyanobacteria sequences. It’s important to filter our these eukaryotic sequences before we start playing with statistical analyses.

After dropping unknown or undesired sequences from our combined ASV data frame, it’s time to create an ASV table that is properly formatted for the vegan package. This ASV table must be a Samples x Species matrix, in which our Sample IDs as our row names and our ASV IDs as our column names.

NOTE: We could have made this ASV table earlier immediately after importing the ASV count data by transposing the table with the t() function. However, I want to have an ASV table that excludes taxa I do not want in my data set, like ASVs attributed to Chloroplast sequences or ASVs attributed to unknown Phyla.

The last thing we need to do before we jump into statistical analyses is reformat the metadata to be in the same order (by rows) as our ASV table. This is a crucial step! Though it may appear minor, certain functions will not correctly analyze your data if your metadata and your ASV table are not arranged in the same order. This next step will only work if the two data frames we are reordering have the same number of rows AND the same row names.

Now that all of our files are in R and correctly formatted, we can start some statistical analyses!

4.2 Alpha Diversity & Species Richness

4.2.1 Rarefaction Curves

To calculate species richness and alpha diversity (using the Shannon-Wiener index), I like to use functions from the vegan package. Before I get to the alpha diversity and species richness, I will calculate a rarefaction curve for my ASV table. The rarefaction curve tells us that after resampling a pool of N individuals per sample (x-axis), we will identify a certain number of species in said sample (y-axis). This can give us an idea if any sample is more/less species rich than other samples, which can be useful to identify outliars.

Figure 7: Rarefaction Curve Example, Colored by Sample Category


In this rarefaction curve, each curve is colored based on its sample category: red represents “Clear Cut Soil”, yellow represents the “Gopher” category , light blue represents the “No Gopher” category, and dark green represents “Old Growth” soil. Based on this rarefaction curve, it appears that samples within the Old Growth category have a smaller, average sample size compared to the other sample categories.

Now that I’ve viewed the rarefaction curve and looked for outliers, we can move onto the alpha diversity and species richness steps.

For alpha diversity and species richness measurements, we are going to use raw data. The use of raw data for any kind of analysis is quite controversial because not all of the samples have the same number of observations: for example, one of our samples may have thousands of ASV counts, whereas other samples can be much smaller or larger than that. Transforming our data can allow us to view the actual distribution of our data, revealing patterns that may have been difficult to observe in the raw data. You can find a helpful example of data transformations and the benefits via its respective Wikipedia page.

Some microbiologists would tell you that we should rarefy our data before moving onto any diversity assessments or downstream analyses. Rarefying is a type of data transformation that involves finding the sample with the minimum number of counts in all of your samples, then scaling all of your sample counts down to this size. As described with the rarefaction curve, rarefying allows you to see 1. the number of species across samples and 2. the abundance of said species across samples when sampling based on a given minimum. Historically, rarefaction was the strategy used to transform microbial data. However, more recently many statisticians have advised against rarefaction as we tend to lose a lot of information regarding low abundance OTUs/ASVs. For more information on why rarefaction is not a useful transformation method, please read “Waste Not, Want Not: Why Rarefying Microbiome Data is Inadmissable” by McMurdie and Holmes (2014).

When we get to the section on calculating beta diversity, I will provide more insight into which transformation method(s) I use and why.

Alpha diversity is a way to measure within-sample diversity, using an equation that considers the richness of certain species as well as the evenness of those species. The vegan package has a diversity() function that allows one to specify which diversity index the user would like to use for calculating alpha diversity. Here I use the Shannon-Wiener index for calculating alpha diversity. In order to calculate Shannon Diversity, we have to calculate Shannon Entropy, then take calculate the exponential value of the Shannon Entropy (e to the power of Shannon Entropy). We can also calculate species richness (i.e., how many species are present in each sample) using the specnumber() function from the vegan package. Once we’ve found species richness and Shannon diversity, we can combine these values into one data frame, then merge this data frame with our metadata to create one dataframe containing: Shannon entropy, Shannon diversity, species richness, and your sample metadata.

We can now use the data frame we made with our alpha diversity, species richness, and our metadata to create some nice figures. First want to ensure that the category of interest (i.e., in this example that will be “group”) is the right class of variable for generating this figure. Because we are using a categorical identifier, it is wise for us to make sure that our Category variable is in the factor format.

Now let’s make some pretty figures with ggplot2! Using ggplot2, we can specify what type of plot to make, the color palette you’ll use, the size(s) of your font, etc. If you’re interested in everything that ggplot2 can do, please check out this amazing ggplot2 Cheat Sheet. We are also using ggpubr, a wrapper for ggplot2 that allows for easy manipulation and export of ggplot figures. For more information on ggpubr, please check out the package website.

Here we are going to create box-and-whisker plots of our alpha diversity and species richness data. The first plot will display the alpha diversity across of our groups, and the second plot will display the species richness of these groups. The y-axis will show the Shannon diversity in the first plot, and the species richness in the second plot. For both plots, the x-axis will display the Group labels.

Each of the individual box-and-whisker plots will be assigned a different color based on the Category variable using the $Category_col variable we created earlier.

Figure 8a: Boxplot of Alpha Diversity by Sample Category


We observe the highest average Shannon diversity within the Clear Cut Soil category, followed by the Gopher category, then the No Gopher category, and the Old Growth category. Though this figure is helpful for comparing these categories to one another, we cannot really glean meaningful statistical inforamtion from this boxplot.

Not only does ggpubr help us with arranging and saving figures, but we can also use some if its functions to add statistics directly onto our figures with ease. In the next code chunk we use a function called stat_compare_means() which allows you to compare different sets of samples to each other. Because we have already assigned our samples to Groups, we can compare the means across our multiple samples.

We can compare the means of each sample to each other in a pair-wise fashion by using either a T-test (t.test) or a using a Wilcoxon test (wilcox.test), or compare the means across all of our samples at once using an Analysis of Variance aka ANOVA (anova) or a Kruskal-Wallis test (kruskal.test). For more information on how to use the stat_compare_means() function to add statistics to your plots, please see this website.

Deciding on whether to use a T-test verses a Wilcoxon test, or an ANOVA verses a Kruskal-Wallis test, depends on whether your data fulfills certain assumptions held by these statistical tests. One of the assumptions for a T-test and an ANOVA is that the data is normally distributed. We can test for normality using the Shapiro-Wilk test. The null hypothesis for the Shapiro-Wilk test is that the data is normally distributed. This means that if your p-value for the Shapiro-Wilk test is > 0.05, then the null hypothesis is accepted and the data is in fact normally distributed. However, if p is < 0.05, then the null hypothesis is rejected and your data are not normally distributed.

We can also use a Q-Q plot to compare our data with a theoretical normal distribution. These plots show the quantiles for our data in the y-axis, and the theoretical quantiles for a normal distribution on the x-axis. If our data points lie on the line in the Q-Q plot, then the data is considered normally distributed. Skewed data will contain points that are further from the line, curving one way or another.

Let’s run a Shapiro-Wilks test using our species richness results, and use a Q-Q plot to see the distribution of these richness values.

NOTE: diversity and richness are usually not normally distributed, but it’s still important to always see how the data are distributed if you plan on running statistical tests that assume normality.

Figure 8b: Normal Q-Q Plot: Species Richness


Because our p-value for the Shapiro-Wilks test is > 0.05, we’ve determined that our species richness values are not normally distributed. Because of this, we will use a Wilcoxon test (rather than a T-test) to compare the means of our sample groups in a pairwise fashion. Because we only have two groups in this example, we cannot run an Kruskal-Wallis test. Kruskal-Wallis tests and ANOVAs are used when comparing three or more groups.

In the boxplot below, I have only included a few pairwise group comparisons as to not overwhelm the plot. If you’d rather use * as indicators of statistical significance instead of using the p-values themselves, you can change the label parameter in the stat_compare_means() function from "p.format" to "p.signif".

Figure 8c: Boxplot of Species Richness by Sample Category


I did not want to overwhelm you with multiple pairwise group comparisons on this figure, so you are only seeing the results of three Wilcoxon tests comparing Clear Cut Soil to Old Growth samples, Gopher to No Gopher samples, and No Gopher to Old Growth Samples. The * indicate significance levels, with * being p <= 0.05, ** being p <= 0.01, *** is p <= 0.001, and **** is p <= 0.0001. The symbol ns stands for not significant. Figure 8c shows us that the Clear Clut Soil samples have a significantly higher average Species Richness than the Old Growth samples. We can also see that our No Gopher samples are significantly higher in species richness than the Old Growth samples, but the difference between the species richness averages in the Gopher verses No Gopher samples is not statistically significant. You can also see a printed value for the Kruskal-Wallis test, which here is p = 0.00008. This test indicates that there is a significant difference in average species richness across all of the sample categories.

4.3 Beta Diversity

Before going any further, we should transform our data. Data transformation helps us to better interpret our data by changing the scales in which we view our data, as well as reducing the impact of skewed data and/or outliers in our data set. We can also perform a transformation that normalizes our data, aka changing the distribution of our data to be a normal (i.e., Gaussian) distribution, which is useful for running certain statistical tests that assume normality (like T-tests, ANOVAs). For more on why you should transform your data and what kind of transformations are out there, check out the resources included in this very helpful Medium article. I have also found this Medium article on log transformations helpful as well.

Two useful transformations I have seen used are the variance stabilizing transformation (i.e, VST) and the centered log-ratio transformation (i.e, CLR). For information on how to employ this particular transformation, please check out this tutorial by the legendary bioinformatician Dr. Mike Lee. Though I won’t be using the VST transformation, I have not found any literature saying that the CLR transformation is better than VST. The CLR transformation appears to be popular among statisticians, which is why I am choosing to go this route.

We will use the microbiome package to CLR transform our count data for creating clustering dendrograms and ordinations. The CLR transformation is recommended in the paper “Microbiome Datasets Are Compositional: And This Is Not Optional” by Gloor et al. (2017), which proposes that microbiome data sets are compositional, meaning they describe relationships between multiple components. Gloor et al. (2017) argues that the reason that CLR transformations are ideal for compositional data is because 1. ratio transformations are useful for detecting relationships in both count data and proportion data, and 2. because log-ratios specifically make the data symmetric in a log space.

A requirement for the microbiome::transform function is that it requires a phyloseq object (i.e., a type of object created by the phyloseq package) as input. Let’s create our phyloseq object with our bac.ASV_counts matrix and add a small pseudocount to the cells in our matrix containing 0s. It’s important to change all of the 0s to a very small number so that the log of 0 is not taken in our transformation. To better understand how CLR transformation works and why its a useful transformation for microbiome data, watch this great YouTube video created by Dr. Thomas Quinn and read “Normalization and microbial differential abundance strategies depend upon data characteristics” by Weiss et al. (2017). For more information on the pros and cons of the CLR transformation, please read “A field guide for the compositional analysis of any-omics data” by Quinn et al. (2019).

Now that we have our transformed ASV counts, we can create a Euclidean distance matrix that will describe how close (aka similar) our samples are to each other based on their microbial composition.

Figure 9: Hierarchical Clustering (with Centered Log-Ratio Transformed Data)


Though there are some samples not clustered within their sample categories, overall it appears that samples from specific categories form distinct clusters. This indicates thats generally, samples from the same category have similar microbial community composition. We can also see that most of the samples in the Clear Cut Soil category are more similar to the Gopher and No Gopher samples compared to the Old Growth samples. Though this dendogram is helpful, it’s not as informative as other visualizations we can do. To learn more about how these sample categories’ microbiomes compare to one another, we can use our Euclidean distance matrix (created from CLR transformed ASV counts) to generate an ordination known as Principal Coordinate Analysis, aka a PCoA.

I am not going to get into the math behind a PCoA, but you can learn more by watching this excellent StatQuest YouTube video and this helpful link that describes what a PCoA is and its uses. If you’re interested in learning more about ordinations in general and the impacts they can have on microbiome data, please read “Uncovering the Horseshoe Effect in Microbial Analyses” by Morton et al. (2017).

Generally a PCoA is used at looking at how similar your samples are to each other, and the variability exhibited by your samples, in a reduced dimensional space. The closer two points are in your ordination, the more similar they are. PCoAs yield multiple axes (i.e., principal components) that capture the variation within your data set and are associated with certain values (i.e., eigenvalues) that represent to the magnitude of the variation for each axis. These eigenvalues are relative representations of how important each axis of variation is for describing the data set.

PCoAs were developed so that we can create these ordinations using distances that are NOT Euclidean, for example like Bray-Curtis dissimilarity distances. PCoAs come from Principal Component Analysis, which is specifically used for Euclidean distances. For more information on PCAs, check out this StatQuest YouTube video as well as this tutorial that compares PCAs to PCoAs.

Let’s generate our PCoA and a PCA and check out the proportion of variance explained by our axes. The reason I am showing you how to generate both is because a PCoA using Euclidean distances is equivalent to a PCA. This is a useful way to confirm that our data is actually Euclidean in nature.

The first axis (PC1) of variation describes 8.94% of the variance in the entire data set. The second axis (PC2) describes 5.69% of the variation. Our PC axes generated by our PCoA are equivalent to our axes generated by the PCA, which is to be expected using Euclidean distances. The first 2-3 axes describe the greatest amount of variation in the data set, and are included in the visualization of the PCoA.

To visualize our PCoA with ggplot2, we have to extract the principal coordinates for each sample across our axes of variation and combine these values with our metadata. Then we can make a PCoA ordination and color each sample ID with our variable of interest (in this case, the Sample Category, aka metadata$Category). We can also include the relative variation for each axis in our x-axis and y-axis labels.

Figure 10: Principal Coordinates Analysis, Colored by Sample Category


From this PCoA we can tell that the microbial community composition in the Old Growth samples are similar to one another, forming a tight cluster in the PCoA (points in dark green). The microbial composition of the Clear Cut Soil samples are also similar to one another (points in red). Interestingly, the Gopher and No Gopher samples cluster together, indicating that their microbial communities are similar to each other, regardless of whether a gopher was introduced to the soil or not. It is important to keep in mind that though we see distinct clusters by sample category, the variation explained by this variable is quite low (Axis 1 - 8.94%, Axis 2 - 5.59%).

So this information is helpful, but we are not sure if our categories are significantly similar/dissimilar from each other. To do this, we first need to check the dispersion (aka variance) of the composition data within each group to see if we can even compare these groups to each other. Basically we are finding the spatial median or the centroid of each group in multivariate space, and calculating the distance from each point to the centroid within a respective group or category. The actual distances are reduced to principal coordinates (as is done in a PCA or PCoA) before the distances from each group are compared. We then can use an ANOVA as well as a Tukey’s Honest Signifcant Difference Test (aka Tukey’s HSD) to statistically compare the group dispersions.

To check out our group dispersions and whether or not they are homogeneous (equal/similar to each other in their variance), we will use the betadisper() function from the vegan package. We can then compare the axes of dispersion with the anova() function and compare the spatial means with the TukeyHSD() function.

The ANOVA results tell us that our dispersions by category are significantly different than each other (p=0.0001394), meaning the variance within each category is not homogenous. We can visualize this comparison as well via an ordination (calculated by betadisper()) and a boxplot based on the distance to the centroid for each group.

Figure 11a: Principal Coordinates Analysis of


Figure 11b: Principal Coordinates Analysis. Colored by Sample Category


The reason that our dispersion results are problematic is that if we try to compare the groups using a Permutational Analysis of Variance (aka a PERMANOVA), the significant differences we may see between groups could be attributed to their unequal variances (i.e., dispersion effects) rather than actual differences in the community compositions by category.

A PERMANOVA is similar to an ANOVA in that both analyses compare differences between groups (using sum-of-squares), but a PERMANOVA runs multiple permutations to compare these distances to each other - whereas an ANOVA is comparing group averages to each other without the use of permutations. Another difference is that while ANOVAs assume that the data is normally distributed, the PERMANOVA assumes that the groups have equal variance (dispersion). For more information on PERMANOVAs and comparing group variances, check out this very helpful link by Dr. Joshua Ebner.

Even though we should not run a PERMANOVA with these data, let’s go over how we would run a PERMANVOA. The most crucial thing about running a PERMANOVA in R is that your feature table and your metadata need to be in the same order by row! The program does not know to match up sample IDs or labels to each other, so you have to confirm that your data frames are arranged in the same way by row before running the PERMANOVA. We can then run our PERMANOVA including multiple variables of interest. We can see if there are interactions between multiple variables and our compositional data respectively with + (i.e, var1 + var2), or we can check for interactions between our variables and our composition data with * (i.e, var1 * var2).

# check your rownames and the order of the rownames for ordering step
rownames(b.clr)
rownames(metadata)
# in case you need to reorder
b.clr=b.clr[rownames(metadata),]
# sanity check
rownames(b.clr)
rownames(metadata)

# PERMANOVA requires assumption of homogenous within-group disperions to ensure observed differences in groups are real
perm1<-adonis2(b.clr ~ Category, data = metadata, permutations = 999, method="euclid", by='terms') # looks for interactions between predictor variables
perm1

perm2<-adonis2(b.clr ~ Category*layer, data = metadata, permutations = 999, method="euclid", by='terms')
# * -> for interactions between predictor variables; + -> interactions with multiple variables but not between them or combined interactions
perm2
adonis2(b.clr ~ Category+layer, data = metadata, permutations = 999, method="euclid", by='terms')

# export PERMANOVA results to csv
perm1_results<-data.frame(DF=perm1$Df, SumofSqs=perm1$SumOfSqs, R2=perm1$R2, F=perm1$F, p=perm1$`Pr(>F)`)
rownames(perm1_results)<-rownames(perm1)
perm1_results
write.csv(perm1_results,"16S_PERMANOVA_Results.csv",row.names=TRUE)

# save PERMANOVA results in a nice table
tab <- ggtexttable(perm1_results, theme = ttheme("light"))
tab2<- tab %>%
  tab_add_title(text = "PERMANOVA Results: Composition ~ Category", face = "bold", padding = unit(1, "line")) %>%
  tab_add_footnote(text = "Using Euclidean Distance of CLR-Transformed Data (Aitchison Distance)", size = 10, face = "italic")

# save table as png
png('permanova_table_test.png',width = 700, height = 600, res=200)
tab2
dev.off()
Figure 12: Table of PERMANOVA Results (by Category)


As predicted by the betadisper() results, we are seeing a significant difference in community composition between our groups. Again, for this data set we can’t really know if these differences are meaningful because the within-group disperions (aka variances) are NOT homogenous - so the significance we are seeing here is likely due to dispersion effects rather than real differences between groups.

4.4 Relative Abundance

So we know our microbial communities are similar to each other if they’re from the same sample category, and that the microbiomes from each category are quite different from each other…but what about the actual taxa within those communities? Calculating the relative abundance allows us to determine how abundant certain taxa are compared to the rest of the taxa within a specific sample. In order to get a better sense of the microbial community in our samples, we need to calculate the relative abundance at varying taxonomic levels. We can also calculate the relative abundance by group; for example, we can determine the relative contribution of taxa 1 compared to the rest of the microbial commiunity in Group A.

Let’s merge our ASV taxonomy table to our metadata table. This will help us create a data frame that we can use to get the sum of our ASV counts across specific taxonomic levels.

Using the all_bac data frame containing our data and metadata altogether and the dcast() function, we can create data frames that contain counts by taxonomic levels. We will then use these data frames to get the relative abundance at specific taxonomic levels. Because I am using a bacterial data set, I am only going to calculate the relative abundance at the phyla and class levels. Though we could get deeper taxonomic resolution, the visualization at these levels can be really overwhelming because there are SO many taxa present.

Remember that when calculating your relative abundance, each count per taxa in each sample is divided by the sum of all the counts for that sample. This means that when you sum up all of the relativized values for each sample, they should sum up to 1. We can use this fact to help us check whether our relativization worked.

Now that we have our relativized phyla counts, let’s visualize the relativized counts by phyla using a stacked barplot.

Figure 13: Relative Abundance by Phyla Across Samples


We can see that across samples, the two most relatively abundant phyla are Abditibacteriota in coral, and Proteobacteria in light purple. Though this is helpful, it’s hard to determine which categories our samples belong to.

Now that we have calculated the relativized counts by class according to our sample categories, let’s visualize it! Instead of a stacked barplot, we could visualize this in a different way - shout out to Dr. Mike Lee for sharing this technique in his amplicon tutorial!

Dr. Lee visualized the proportions of gene copies recovered in a way that he calls taxonomic summaries. These summaries allow us to view the relative abundance (or gene copy number, or any abundance measurement) of our microbes by having our x-axis be the taxa themselves rather than the sample IDs or the group labels. Personally I find these taxonomic summaries useful for comparing the relative abundance of specific taxa across multiple groups. To provide an example that’s easy to look at, we are going to view microbial classes by category that have a relative abundance of at least 1% or higher per sample.

Figure 14: Taxonomic Summary of Bacterial Relative Abundance (> 1%) by Class Per Sample. Colored by Sample Category.


This taxonomic summary tells us that Old Growth samples have a higher relative abundance (>20-40%) of the Acidimicrobiia class compared to samples belonging to the other categories. Gammaproteobacteria appears to be the most abundant class in Clear Cut Soil samples, with some samples having a relative abundance of this taxa that is greater than 20%. Interestingly, these samples from the Clear Cut Soil category are the only samples contain Bacilli, though this could be due to our 1% relative abundancer threshold for this figure.

Maybe it would be helpful to know the relative abundance of taxa by the categories themselves rather than the relative abundance by sample. We can easily calculate this and visualize a stacked barplot to compare sample categories.

head(all_bac)
# by class + elevation
bac.cls.cat <- as.data.frame(dcast(all_bac,Category~Class, value.var="Counts", fun.aggregate=sum)) ### 
head(bac.cls.cat) # counts by class + elevation
rownames(bac.cls.cat)<-bac.cls.cat$Category

b.RA_cls.cat<-data.frame(decostand(bac.cls.cat[,-1], method="total", MARGIN=1, na.rm=TRUE))  
# relative abundance of taxa data where everything is divided by margin total (default MARGIN = 1 = rows) -- rows = samples
rowSums(b.RA_cls.cat) # sanity check
b.RA_cls.cat$Category<-rownames(b.RA_cls.cat)
head(b.RA_cls.cat)

#melt down relativized data to merge with metadata
b.cls.cat_m<-melt(b.RA_cls.cat, by="Category")

head(b.cls.cat_m)
colnames(b.cls.cat_m)[which(names(b.cls.cat_m) == "variable")] <- "Class"
colnames(b.cls.cat_m)[which(names(b.cls.cat_m) == "value")] <- "Counts"
head(b.cls.cat_m) ## relative abundance based on sum of counts by class!

# Subset data to include only classes with relative abundance > 1%
b.cls.cat.1<-subset(b.cls.cat_m, b.cls.cat_m$Counts>0.01)

c1<-ggplot(b.cls.cat.1, aes(x=Category, y=Counts, fill=Class))+geom_bar(stat="identity",colour="black")+scale_x_discrete(labels=c("ClearCutSoil"="Clear Cut Soil", "Gopher"="Gopher", "NoGopher"="No Gopher", "OldGrowth"="Old Growth"))+theme_classic()+
  labs(title = "Microbial Class Relative Abundance", x="Sample Category", y="Relative Abundance", fill="Class", subtitle="Only Including Taxa with Relative Abudance > 1%")+
  theme(axis.title.x = element_text(size=13),axis.title.y = element_text(size=13),axis.text = element_text(size=11),axis.text.x = element_text(vjust=1),legend.title.align=0.5, legend.title = element_text(size=13),legend.text = element_text(size=11),plot.title = element_text(size=15),plot.subtitle = element_text(size=13))

ggsave(c1,filename = "16S_Class_1percent_RA_byCategory.png", width=12, height=10, dpi=600)
Figure 15: Relative Abundance (> 1%) of Bacterial Classes by Sample Category


Instead of viewing the relative abundance of classes by sample as we did in the taxonomic summary, we are now looking at the relativea abundance of taxa within the categories themselves. For example, we can tell that out of all of the categories, samples within the Old Growth category have the highest relative abundance of the Acidobacteriae class. Gammaproteobacteria appears to be the most abundant in the Clear Cut Soil category, whereas Bacteroidia is the most abundant class in the No Gopher category. Though our taxonomic summaries helped us to see these trends, these stacked bar plots offer some more insight to category-specific trends.

4.5 Environmental Drivers

The last step in this workflow is to determine how our microbial community composition corresponds to environmental variables. To do this, we can run either a Canonical Correspondence Analysis (aka CCA) or a Redundancy Analysis (aka RDA). These are ordination techniques that allow us to correlate our composition data with our environmental variables to explain the variance in our microbial community structures. CCAs and RDAs utilize eigenanalysis, which is quite out of the scope of this tutorial. I highly recommend watching this YouTube video detailing what eigenvectors and eigenvalues are and why they’re essential in linear algebra and statistics.

A CCA considers how community composition relates to a set of defined constraints, which would be your environmental variables of interests. CCAs assume that your composition data will have a unimodal (aka have one maximum value) relationship with the environmental variables. Typically for a CCA you already know exactly what environmental variables you will be considering. It’s important to note that if your model (how your environmental variables influence your community composition) is not statistically significant, then your results are not meaningful enough to discuss.

A RDA can tell us whether our environmental variables explain the variation we observe in our microbial communities, or whether or not their impact on community structure is redundant (i.e., the variation is explained by something else). RDAs assume that your composition data will have a linear relationship with your environmental data. Variance described by our environmental response variables is known as constrained variance, whereas variance described by unknown variables (or variables excluded from our response variables) is called unconstrained variance.

So at this point you’re probably asking yourself “how do I know if my data has a unimodal or linear relationship with my environmental variables?”, and that’s a great question! We can use a Detrended Correspondence Analysis (aka DCA) to help us decide whether we should use a CCA or an RDA for our composition data. If the length of our first DCA axis is longer than 4, that indicates that our data is heterogeneous and requires unimodal methods like a CCA. If the length of our first DCA axis is less than 3, the data is considered homogeneous and thus requires linear methods like an RDA. If your first axis length is between 3 and 4, then it’s up to you to choose whether to use an RDA or a CCA.

For more information on how CCAs and RDAs are calculated and their applications, please view this helpful tutorial for a deep look into different Euclidean-based ordinations. I also recommend watching Dr. Ralf Shaefer’s YouTube video on RDAs as well as this helpful YouTube video that walks you through different ordination techniques using the vegan package in R.

Let’s first run a DCA to determine the length of our first axis and decide which ordination technique to use. For this we are going to use our CLR-transformed ASV table as our site x species matrix. The decorana() function from vegan performs the DCA, which requires that the row sums are greater than zero. Because clr transformed data can have negative values, I’ve decided to add a pseudocount of 1 to the matrix to fulfill the requirements of the decorana() function.

Note: The chemical data I will be using for this portion of the workflow is not real data from this project. These data come from a different project and are being used just for this example.

Our first DCA axis has a length of 1.7990, which suggests that we should run an RDA with our data.

We first want to scale our environmental data so that they are comparable to each other, despite their initial measurements being vastly different (i.e., your mg of Cu and mg of P are not immediatley comparable). At this point, we could jump into doing an RDA…but it may be more informative if we look at the environmental drivers of each category or site one by one rather than considering the drivers across sites.

To subset our data, we will use the lapply() function: if metdata$Category is the same as (==) the element(s) in site_list, subset the metadata data frame by the site_list element into a new list called site_subsets. After this step, we will have a list where each element in the list (site_subsets[[i]]) will contain all of the metadata for each respective Category.

I have also provided examples on how to subset lists for practice - indexing in R can be tricky in the beginning, and personally I need to review it often.

The next step is to take each element in our list and make a new data frame, so that we have data frames of the metadata separated by Category. I have written a custom function called df_specific.subset to do this, but if you know an easier way of doing this, feel free to use that instead!

At the end of this step, we should have data frames of each category containing all of their respective metadata.

Now that we have data frames containing all of the metadata by category, we need to match these data frames up with our composition data (our transformed ASV table). We are going to use another custom function I wrote called match_dat to do this. First we have to run the function to make sure it’s in our Global Environment, and then we can use it. To save time, I am only going to use one of our category data frames to match up with our transformed composition data.

Now that we have the metadata and the corresponding transformed ASV counts for the samples in the Clear Cut Soil category, we can run an RDA with the rda() function from the vegan package to see if our variables of interest drive our microbial community composition.

If your DCA 1 axis was longer than 3 and you needed to run a CCA, all of these steps would be the same except that instead of using the rda() function by vegan, you would use the cca() function. Similar to the decorana function, you may need to add a small pseudocount to your transformed feature table before you run the cca() function.

To view the results of our rda, we will use the summary() function, and we can check the amount of variation explained by the model (i.e., R2) using the RsquareAdj() function. It’s wise to use the adjusted R2 value because the more variables included in the model, the more inflated your R2 value will be. We can also check out the individual terms of the model and their respective impact via an anova(). Lastly, we can see exactly which variables are the drivers of variation in our model using the ordistep() function. Based on the results of this last step, we can adjust our model and re-run the RDA.

Because our model was not significant, normally we would stop here and not visualize this RDA. However, the point of this workflow is to provide you with the code to run and visualize these analyses yourselves, so I am going to show you how I would visualize this RDA.

The ggvegan package is a useful package that takes advantage of some of the difficulties of ggplot2 when it comes to plotting ordinations. This package has a function called autoplot which helps us easily plot an RDA.

Figure 16: Redundancy Analysis Autoplot by ggvegan


We can also extract data from the summary(rda) object and use this information to create our own RDA plot with ggplot2. I am going to show you two versions of the plot: one version that has not been altered in any way (Figure 17a), and an RDA where the axes pointing to our environmental variables have been amplified to match the autoplot (Figure 17b).

Figure 17a: Redundancy Analysis w/ ggplot2 Example 1


Figure 17b: Redundancy Analysis w/ ggplot2 Example 2


The first ggplot does not amplify the biplot axes from our RDA as we see with the RDA created by autoplot(). The second plot looks more like what autplot() produced, which contains axes that were amplified to clearly show the trends. It’s important to be cautious of this type of ordination manipulation because you do not want to “identify” relationships in your data that do not actually exist. That is why I recommend comparing your RDA axes from the RDA summary() (i.e., your site, biplot, and species axes) so that you know exactly how much variation is being described by your environmental variables.

If our RDA had shown these variables as being significant drivers of community composition, that is what I would describe here. It is also wise to include a table that describes the significance of your RDA or CCA as well as the variation (i.e., adjust R2 value) explained by your model. However, because the chemical data used for this portion of the workflow is fictional, I won’t describe this figure in depth.

Because none of my environmental variables were significant drivers of community composition, I should not have included them in my RDA. Remember that it’s important to only show environmental variables on your RDAs or CCAs that significantly describe the variation in your compositional data that’s explained by your model.

5 Conclusion

Congratulations, we have reached the end of this workflow! By now you should have some analyses and figures for your data set, and you can start to get to work on interpreting your results. We can do ourselves a favor and save everything in our global environment by running the following line of code.

Please keep in mind that this is a collection of scripts and analyses that I have created, and this is by no means a Bible for microbial ecology! You do not have to take all of these steps with your own data, and if you have more efficient ways of doing some of these steps, then I encourage you to use those methods. This worfklow is meant to serve as a jumping off point rather than a final destination, so please do whatever you think is necessary to investigate the patterns you’re observing in your microbial communities.

If you’d like to get a hold of me to offer me feedback about this workflow, or if you’d like to discuss more of these steps/programs/stats in depth, do not hestitate to reach out. My contact information is in the About Me section.

Thank you for following along!

6 Version Information

## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## Random number generation:
##  RNG:     Mersenne-Twister 
##  Normal:  Inversion 
##  Sample:  Rounding 
##  
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] devtools_2.4.2           usethis_2.0.1            microbiome_1.12.0        shades_1.4.0            
##  [5] nationalparkcolors_0.1.0 wesanderson_0.3.6        reshape2_1.4.4           reshape_0.8.8           
##  [9] tidyr_1.1.3              readxl_1.3.1             dplyr_1.0.7              plyr_1.8.6              
## [13] ape_5.5                  scales_1.1.1             ggpubr_0.4.0             vegan_2.5-7             
## [17] lattice_0.20-44          permute_0.9-5            ggplot2_3.3.5            phyloseq_1.34.0         
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15          colorspace_2.0-2    ggsignif_0.6.3      ellipsis_0.3.2      rio_0.5.27         
##   [6] rprojroot_2.0.2     XVector_0.30.0      fs_1.5.0            rstudioapi_0.13     remotes_2.4.0      
##  [11] fansi_0.5.0         codetools_0.2-18    splines_4.0.5       cachem_1.0.6        knitr_1.34         
##  [16] pkgload_1.2.2       ade4_1.7-17         jsonlite_1.7.2      broom_0.7.9         cluster_2.1.2      
##  [21] compiler_4.0.5      backports_1.2.1     assertthat_0.2.1    Matrix_1.3-4        fastmap_1.1.0      
##  [26] cli_3.0.1           htmltools_0.5.2     prettyunits_1.1.1   tools_4.0.5         igraph_1.2.6       
##  [31] gtable_0.3.0        glue_1.4.2          Rcpp_1.0.7          carData_3.0-4       Biobase_2.50.0     
##  [36] jquerylib_0.1.4     cellranger_1.1.0    vctrs_0.3.8         Biostrings_2.58.0   rhdf5filters_1.2.1 
##  [41] multtest_2.46.0     nlme_3.1-153        iterators_1.0.13    xfun_0.26           stringr_1.4.0      
##  [46] ps_1.6.0            openxlsx_4.2.4      testthat_3.0.4      lifecycle_1.0.0     rstatix_0.7.0      
##  [51] zlibbioc_1.36.0     MASS_7.3-54         hms_1.1.0           parallel_4.0.5      biomformat_1.18.0  
##  [56] rhdf5_2.34.0        yaml_2.2.1          curl_4.3.2          memoise_2.0.0       gridExtra_2.3      
##  [61] sass_0.4.0          stringi_1.7.4       S4Vectors_0.28.1    desc_1.3.0          foreach_1.5.1      
##  [66] BiocGenerics_0.36.1 pkgbuild_1.2.0      zip_2.2.0           rlang_0.4.11        pkgconfig_2.0.3    
##  [71] evaluate_0.14       purrr_0.3.4         Rhdf5lib_1.12.1     cowplot_1.1.1       tidyselect_1.1.1   
##  [76] processx_3.5.2      magrittr_2.0.1      R6_2.5.1            IRanges_2.24.1      generics_0.1.0     
##  [81] DBI_1.1.1           pillar_1.6.2        haven_2.4.3         foreign_0.8-81      withr_2.4.2        
##  [86] mgcv_1.8-36         survival_3.2-13     abind_1.4-5         tibble_3.1.4        crayon_1.4.1       
##  [91] car_3.0-11          utf8_1.2.2          rmarkdown_2.11      data.table_1.14.0   callr_3.7.0        
##  [96] forcats_0.5.1       digest_0.6.27       stats4_4.0.5        munsell_0.5.0       bslib_0.3.0        
## [101] sessioninfo_1.1.1

7 About Me

My name is Hannah and my pronouns are they/them. I am currently a PhD Student at UC Riverside in the Genetics, Genomics, and Bioinformatics PhD program and a member of Dr. Emma Aronson’s lab.

If you have any questions regarding this workflow and the scripts I used, do not hesitate to contact me. Or, if you’d like to talk bioinformatics and all things ’omics, I would love that too!

If you’d like to learn more about my research background and check out some of my other work, feel free to visit my website. It is a work in progress, just fyi!

For some personal information about me, check out my page with 500 Queer Scientists, and my twitter if you’re curious.

Thank you so much for checking out my workflow!
website

References

Andrews, S. n.d. “FastQC: A Quality Control Tool for High Throughput Sequence Data.”

Bukin, Yu S., Yu P. Galachyants, I. V. Morozov, S. V. Bukin, A. S. Zakharenko, and T. I. Zemskaya. 2019. “The effect of 16s rRNA region choice on bacterial community metabarcoding results.” Scientific Data 6. The Author(s): 1–14. https://doi.org/10.1038/sdata.2019.7.

Bushnell, Brian. n.d. “BBMap.” sourceforge.net/projects/bbmap/.

Callahan, Benjamin J., Paul J. McMurdie, and Susan P. Holmes. 2017. “Exact sequence variants should replace operational taxonomic units in marker-gene data analysis.” ISME Journal 11 (12). Nature Publishing Group: 2639–43. https://doi.org/10.1038/ismej.2017.119.

Callahan, Benjamin J., Paul J. McMurdie, Michael J. Rosen, Andrew W. Han, Amy Jo A. Johnson, and Susan P. Holmes. 2016. “DADA2: High-resolution sample inference from Illumina amplicon data.” Nature Methods 13 (7). Nature Publishing Group: 581–83. https://doi.org/10.1038/nmeth.3869.

Edgar, Robert C., and Henrik Flyvbjerg. 2015. “Error filtering, pair assembly and error correction for next-generation sequencing reads.” Bioinformatics 31 (21): 3476–82. https://doi.org/10.1093/bioinformatics/btv401.

Gloor, Gregory B., Jean M. Macklaim, Vera Pawlowsky-Glahn, and Juan J. Egozcue. 2017. “Microbiome datasets are compositional: And this is not optional.” Frontiers in Microbiology 8 (NOV): 1–6. https://doi.org/10.3389/fmicb.2017.02224.

Illumina. n.d. “High-Speed, Multiplexed 16S Microbial Sequencing on the MiSeq System.” http://files/998/High-Speed, Multiplexed 16S Microbial Sequencing o.pdf.

Lee, Michael. 2019. “Happy Belly Bioinformatics: an open-source resource dedicated to helping biologists utilize bioinformatics.” Journal of Open Source Education 2 (19): 53. https://doi.org/10.21105/jose.00053.

McMurdie, Paul J., and Susan Holmes. 2014. “Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible.” PLoS Computational Biology 10 (4). https://doi.org/10.1371/journal.pcbi.1003531.

Morton, James T., Liam Toran, Anna Edlund, Jessica L. Metcalf, Christian Lauber, and Rob Knight. 2017. “Uncovering the Horseshoe Effect in Microbial Analyses.” mSystems 2 (1): 1–8. https://doi.org/10.1128/msystems.00166-16.

Prodan, Andrei, Valentina Tremaroli, Harald Brolin, Aeilko H. Zwinderman, Max Nieuwdorp, and Evgeni Levin. 2020. “Comparing bioinformatic pipelines for microbial 16S rRNA amplicon sequencing.” PLoS ONE 15 (1): 1–19. https://doi.org/10.1371/journal.pone.0227434.

Quinn, Thomas P., Ionas Erb, Greg Gloor, Cedric Notredame, Mark F. Richardson, and Tamsyn M. Crowley. 2019. “A field guide for the compositional analysis of any-omics data.” GigaScience 8 (9). Oxford University Press: 1–14. https://doi.org/10.1093/gigascience/giz107.

Rausch, Philipp, Malte Rühlemann, Britt M. Hermes, Shauni Doms, Tal Dagan, Katja Dierking, Hanna Domin, et al. 2019. “Comparative analysis of amplicon and metagenomic sequencing methods reveals key features in the evolution of animal metaorganisms.” Microbiome 7 (1). Microbiome: 1–19. https://doi.org/10.1186/s40168-019-0743-1.

Vargas-Albores, Francisco, Luis Enrique Ortiz-Suárez, Enrique Villalpando-Canchola, and Marcel Martínez-Porchas. 2017. “Size-variable zone in V3 region of 16S rRNA.” RNA Biology 14 (11): 1514–21. https://doi.org/10.1080/15476286.2017.1317912.

Wang, Qiong, George M. Garrity, James M. Tiedje, and James R. Cole. 2007. “Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy.” Applied and Environmental Microbiology 73 (16): 5261–7. https://doi.org/10.1128/AEM.00062-07.

Weiss, Sophie, Zhenjiang Zech Xu, Shyamal Peddada, Amnon Amir, Kyle Bittinger, Antonio Gonzalez, Catherine Lozupone, et al. 2017. “Normalization and microbial differential abundance strategies depend upon data characteristics.” Microbiome 5 (1). Microbiome: 1–18. https://doi.org/10.1186/s40168-017-0237-y.